%%% Estimating Number of Common Factors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Bai and Ng (2002) criteria applied to errors from
%%%% Prewhitened and Projected Principal Components
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function m = ppbn(x,mm,opt,prop)

n = size(x,2);
if prop==1; xx = prwp(x); %homogenous AR(1) prewhitenning proc.
elseif prop==0; xx = x;
end

m = pcms(xx,mm,opt);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Procedure for estimating number of common factors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [mf] = pcms(x,km,opt)

[t,n] = size(x);

msv = zeros(km+1,6);  

if n < t; [b,c,em,vnt] = pcn(x,km);
    for ki = 0:km
        [f,l,e,vnt] = pcn(x,ki);
        [msv(ki+1,1) msv(ki+1,2) msv(ki+1,3) msv(ki+1,4)...
        msv(ki+1,5) msv(ki+1,6),g,h,j] = bnms(e,ki,em,km);
    end
    
else [b,c,em,vnt] = pct(x,km);
    for ki = 0:km
        [f,l,e,vnt] = pct(x,ki);
        [msv(ki+1,1) msv(ki+1,2) msv(ki+1,3) msv(ki+1,4)...
        msv(ki+1,5) msv(ki+1,6),g,h,j] = bnms(e,ki,em,km);
    end

end

if opt==2;
    mf = ceil(mean(minindc(msv(:,5))-1));
elseif opt==1;
    mf = ceil(mean(minindc(msv(:,1:2))-1));
elseif opt==3;
    mf = ceil(mean(minindc([msv(:,1:2) msv(:,4:5)])-1));
end

if mf == 0; mf = 1; end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Bai and Ng (2002) factor model selection criteria
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [pc1, pc2, pc3, ic1, ic2, ic3, bic1, bic2, bic3] = bnms(e,k,ek,km)

[t,n] = size(e);
w = sum(sum(e.*e))/(n*t);
c = min(sqrt(n),sqrt(t));

wk = sum(sum(ek.*ek))/(n*t);

pc1 = w + wk*k*((n+t)/(n*t))*log((n*t)/(n+t));
pc2 = w + wk*k*((n+t)/(n*t))*log(c^2);
pc3 = w + wk*k*(log(c^2)/c^2);
ic1 = log(w) + k*((n+t)/(n*t))*log((n*t)/(n+t));
ic2 = log(w) + k*((n+t)/(n*t))*log(c^2);
ic3 = log(w) + k*(log(c^2)/c^2);
bic1 = w + wk*k*(log(t)/t);
bic2 = w + wk*k*(log(n)/n);
bic3 = w + wk*k*((n+t-k)*log(n*t))/(n*t);
